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Abstract 


At high Mach numbers the two-dimensional time-developing mixing layer develops 
shock waves, positioned around large-scale vortical structures. A suitable numerical 
method has to be able to capture the inherent instability of the flow, leading to the 
roll-up of vortices, and also must be able to capture shock waves when they develop. 
Standard schemes for low speed turbulent flows, for example spectral methods, rely on 
resolution of all flow-features and cannot handle shock waves, which become too thin 
at any realistic Reynolds number. The objective of this work is to study the perfor- 
mance of a class of second-order explicit total variation diminishing (TVD) schemes 
on a compressible mixing layer problem. The basic idea is to capture the physics of 
the flow correctly, by resolving down to the smallest turbulent length scales, without 
resorting to turbulence or sub-grid scale modeling, and at the same time capture shock 
waves without spurious oscillations. The present study indicates that TVD schemes can 
capture the shocks accurately when they form, but without resorting to a finer grid have 
poor accuracy in computing the vortex growth. The solution accuracy also depends on 
the choice of limiter. However a large number of grid points are in general required 
to resolve the correct vortex growth. This phenomenon in computing time-dependent 
problems containing shock waves as well as vortical structures is partly due to the in- 
herent shock-capturing property of all TVD schemes. In order to capture shock waves 
without spurious oscillations these schemes reduce to first-order near extrema and indi- 
rectly produce ‘clipping phenomena’, leading to inaccuracy in the computation of vortex 
growth. Accurate simulation of unsteady turbulent fluid flows with shock waves will 
require further development of efficient, uniformly higher than second-order accurate, 
shock-capturing methods. 
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1. Introduction 


Interest in supersonic combustion has led to renewed efforts to understand the 
physics of high-speed turbulent flows. Mixing in supersonic combustors occurs in com- 
pressible free shear layers, where there can be shock waves and large scale vortical struc- 
tures, which arise from instabilities in the flow. Direct numerical simulations of such 
flows will greatly aid in understanding compressible turbulence, and modeling mixing 
and chemical reactions in high speed flows. The two-dimensional time- developing mix- 
ing layer is a simple prototype of a free shear layer and has been selected for the current 
study, where we investigate the merits of a class of TVD shock-capturing methods. 

The time-developing mixing layer, shown schematically on figure 1, consists of two 
streams of fluid moving in opposite directions, with a smooth velocity profile in between. 
The mean velocity profile has an inflection point and is unstable to small disturbances. 
The linear stability characteristics of the flow have been presented by Sandham and 
Reynolds (1989). It was found that compressibility damps the two-dimensional growth 
of disturbances, and that three-dimensional waves become important at high Mach 
numbers. The physics of the flow are described in detail in Sandham (1989), where it 
is shown that there is a range of free-stream Mach numbers around 0.8 where the two- 
dimensional disturbance is still amplified, leading to a large-scale vortical structure. 
Flow around this vortex is accelerated to locally supersonic speeds and, in order to 
reach the stagnation point in between the vortices, has to compress through a nearly- 
normal weak ( M % 1.2) shock wave. In this flow regime there is an instability to 
be captured, and shock waves to be handled, which represents a severe test for any 
numerical method. The approach in Sandham (1989) was to use spectral and high- 
order Pade finite differences (Lele, 1989) and rely on resolution of all relevant length 
scales. In order to capture shock waves this meant that low Reynolds numbers had to be 
used so that, with 7-8 grid points within the shock wave, there are no oscillations around 
the shock - in fact the computed result is hopefully a full solution of the compressible 
Navier-Stokes equations, including shock wave structure. Those computations were 
limited to a small range of Reynolds numbers - high enough so that the instability 
rolled up strongly enough to form shock waves, but low enough so that the shock waves 
could be fully resolved. This becomes prohibitively expensive as the Reynolds number 
is increased, since it is very inefficient to have to resolve the interned shock structure 
in all the computations. The smallest length scale in the shock wave is of the order of 
the molecular mean free path, while the smallest turbulent length scale is much bigger 
than this. For the transitioned mixing layer the smallest length scale to be resolved is 
a Taylor diffusion scale, where strain balances diffusion, while for a general turbulent 
flow the smallest scale would be of the order of the Kolmogorov scale. 

Computations of the compressible mixing layer using TVD schemes were presented 
by Soestrisno et al. (1988), where the inviscid equations were solved. These compu- 
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tations showed the correct basic structure, but resolution was low. The spatially- 
developing mixing layer was attempted by Guirguis et ai (1987), using an inviscid 
flux-corrected transport (FCT) method. However the computations were at a very high 
convective Mach number where the instability is very weakly amplified. The resulting 
flow fields show low resolution of flow features, indicating that an insufficient number 
of grid points were used in the computations. For turbulent flow computation, and 
calculation of mixing, we need to resolve properly the diffusion length scales in the flow. 

In this paper we assess the performance of a class of second-order explicit TVD 
methods for the viscous time-developing mixing layer, and identify the advantages and 
disadvantages of these methods. The basic idea is to try to resolve down to the smallest 
turbulent length scales, and so capture the physics of the flow correctly, without resort 
to turbulence or sub-grid scale modeling, and at the same time capture shock waves 
without having spurious oscillations. 


2. Problem Formulation 


The governing equations are the two-dimensional time-dependent Navier-Stokes 
equations, with no turbulence model. 

dU dF dG_ J_ 

di + dx dy Re 

where U is (p,pu,pv, e) T . The velocity components in the x and y directions are u and 
v respectively. The density is p and the total energy per unit volume is e. F and G are 
the Euler fluxes in the i and y directions, and F v and G v are the fluxes due to viscous 
terms. Re is the Reynolds number of the flow. The complete formulation in generalized 
coordinates is the same as described in Yee (1987). The mean profiles of velocity and 
temperature are specified by the following relations: 


dF v dG „ 
+ 


dx dy 


(1) 


u = 0.5tanh(2y) (2) 

T = l + M 2 ^-^(1 -u 2 ). ( 3 ) 

The temperature profile satisfies the compressible boundary-layer energy equation, as- 
suming a Prandtl number of unity. The free-stream Mach number is M , and in the 
current work we assume that the free- streams have equal and opposite velocities and 
equal temperatures. Pressure is assumed uniform initially so the mean density profile 
can be obtained directly from the temperature profile. The viscosity p, is assumed to 
follow Sutherlands Law: 


(C 2 /Cl) 1 ' 5 (l + 110.3/Tref) 

“ c 2 /c] + iio.3/r r , f 


(4) 
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where c is the sound speed, subscript 1 refers to the upper free-stream and the reference 
temperature T re f is set to 300 K. The various non-dimensional parameters in the flow are 
the free-stream Mach number, Reynolds number, Prandtl number and Schmidt number. 
The Reynolds number is based on the initial vorticity thickness of the mixing layer, the 
velocity difference across the layer, and the free-stream viscosity. The Prandtl number 
is set to unity, as is the Schmidt number which is required for calculation of diffusion 
terms in the passive scalar equation that is solved alongside the continuity, momentum 
and energy equations. The passive scalar is initialized with a hyperbolic tangent profile 
and tags fluid from each of the free-streams separately, so that mixing of two species in 
the layer can be visualized. 

The flow is periodic in the streamwise direction, x, so periodic boundary conditions 
are applied. Slip- wall boundary conditions are imposed in the normal direction, y , where 
a stretched mesh maps the walls to a large enough distance from the centerline that 
they do not interfere with the computation. 

The flow is unstable to small disturbances. This instability is forced, by adding 
small perturbations, u' and v', to the initial mean flow. 

u' = amp^- sin(27rx/£)e -v ^ 10 (5) 

107T 

v = ampcos(27ra:/£)e -va / 10 (6) 

where ‘amp’ is the amplitude of the disturbance. The perturbation shape on the normal 
velocity component in the y direction is made to resemble the linear eigenfunction, and 
decays to zero in the free-stream. The perturbation in the streamwise velocity is chosen 
so that the entire disturbance is divergence-free. 

In this paper, the disturbance amplitude is chosen to be 0.05 and the wavelength 
L is set to 20 initial vorticity thicknesses. The Mach number is 0.8 and the Reynolds 
number 400. 


3. Numerical Methods 


The Navier-Stokes equations are discretized in conservation law form. The nu- 
merical methods are given here only for one space direction, since terms for the other 
direction are identical. The viscous terms are treated by second-order central differ- 
ences, and for simplicity of presentation are not included here. Let U J be the solution 
at grid point j and at time step n, and let At be the time step and Ax the grid spacing. 
A base method, used for comparison purposes, is the explicit MacCormack method: 

= t/J 1 - A(F” - F/.j) (7a) 
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C r«> = I[^ 1) + p»-A(^ ) 1 -Fj ,) ] (76) 

where X = At/ Ax. Equation (7) can be modified to have the TVD-type property in a 
manner described in Yee (1987,1989). These predictor-corrector TVD-type methods for 
the nonlinear hyperbolic conservation laws (Euler equations) using Roe’s approximate 
Riemann solver can be expressed as: 


rj ,) = U” - X (Fj 1 - Fp.j) 

(8a) 

v f ’ = \ [pj 0 +V7- *(*&>, - Fj”)] 

ST 

00 


(8c) 

When the TVD part (the second term in equation (8c)) is evaluated at the scheme 

is referred to as PCI and when it is evaluated at {7j* it is referred to as PC2. The PC2 
method is proven TVD for the scalar constant coefficient case but requires more storage 
than PCI. See LeVeque and Yee (1988) and Yee (1989) for more details on the behavior 
of PCI and PC2. The matrix R is the right -eigenvector matrix of the flux Jacobian 
A = dFjdU . The specific form used in the current study is shown in Yee (1985). The 
elements of Rj + i are determined by Roe’s averaging (an approximate Riemann solver), 
which can be written as follows: 

Uj + Dua+ i 

u i+\~ 1+D ’ 

(9a) 

Vj + Dv J+ y 

1 + D ’ 

(96) 

Hj + DHj+i 

H > +1 > = 1 + D ’ 

(9c) 


(9d) 

where D = /Eitl and H = - — — . 

V Pi p 

(9e) 

For more details of the development of the schemes the reader is referred to Yee 
(1987,1989). For a symmetric second-order TVD scheme the elements of $ + i, denoted 
b y (*} +f ) s , l = are given by 

(«' +i ) S =K +! l(l-A|«; +s l)(a' +i -Q' +i ) 

(10) 
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where 


a i+\ ~ Rj+i(Uj+i ~ Uj)- 


The value af +1 is the local characteristic speed a 1 evaluated at a symmetric average 
of Uj and Uj+ 1 , using equation (9). For the Euler terms considered here the local 
characteristics are (u — c,u,u,u + c). The following symmetric limiters + i were 

tested: 


Q‘j + J = minmod ( Q !- minmod(a' + } , a( + j ) - a ‘ J+ } , 
Q' + i =minmod(a‘_ i ,a' +} ,a' +5 ), 

(§' +J =mi n mod[ 2 a'. J , 2 a; +§ , 2 a' +| +a' +J )], 


( 12 a) 

(126) 

(12c) 


where the minmod of two arguments x and y is defined as 

minmod(x, y) = sgn(x)max{0,min(|x|,y sgn(x))} (13) 

with sgn(x)=sign(x). These limiters, equations ( 12 a-c), will be referred to as Si, S 2 
and S3 respectively in the test calculations. 

For a second-order upwind TVD scheme the elements of i , denoted by {<j> 1 . A ) u , 
are 3 

(<A' + i) = —<T ( a j+i )(?j+i + 9 j) + (l«J+i + 7 J- 4 I - A ( a 5 + i) 2 ]a5+^ (14a) 


where 


*( z ) = 2^1 ~ 2 ) 


i l i+\ -^(«} + jKf} + , 


(146) 


(14c) 


except when a ( +1 = 0, in which case 7 * A is set to zero. The following upwind limiter 
functions were tested: 

9j = minmod(aJ_i,aj. + i) (15a) 

9 ) ~ ( a j+i a !-i + K+i a i-iO/( a }+i + a )- 1 ) (156) 

9) = i a j- i[( a !+i) + 6 ] + a 5+i[( a !-i ) 2 + 6 ]}/[(aJ + i ) 2 + (a '_ A ) 2 + 28} (15c) 
9) = S • max[0,min(2|aJ + i|,5a'_i),min(|aJ +A |,25aJ._i)];5 = sgn(aJ. + A ) (15d) 

g l j = minmod ( 2 a(_i, 2 a( + A , ^(a( +A +a*_i)). (15c) 
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In the third limiter the constant 6 was chosen to be 10 7 . 

One can modify the limiter easily to obtain a uniformly second-order non-oscillatory 
scheme, which can maintain second-order accuracy at extrema. A variation of the total 
variation bounded (TVB) scheme introduced by Shu (1987) can be obtained (Yee, 1989) 
by redefining the g j function as 

g\ = ^minmod(a' + i,u>aJ_i + MAx 2 sgn(a J+ i )) 

+ ^minmod(a'_i,a?a' + i + MA* 2 sgn(a J _i )). (16) 

In the current implementation u> was set to 1, and the value of M was chosen to be 50, 
following the recommendation of Shu from tests using Burgers’ equation. No attempt 
was made to optimize the values of u and M. 

The above upwind TVD limiters, equations (15a-e), will be referred to as U1-U5, 
and the TVB limiter, equation (16), as U6 in the following sections. Equations (10) 
and (14) are written in the non-entropy satisfying forms since unsteady computations 
are involved. These equations can be easily modified to satisfy an entropy inequality 
(see Yee (1989) for an example). 


4. Results 


The computations were run with a variable time step, at a CFL number, based on 
the convective and viscous terms (MacCormack, 1985), of 0.8. The basic grid used for 
calculation of a single structure was 75 x 75, when 800 time steps were computed. A 
uniformly spaced grid was used in the x direction with box length 20, and a geometric 
stretching was used in the y direction to map the walls to about y = ±50. For com- 
parison purposes we first present results for the growth of a single vortical structure at 
M = 0.8 on a 75 x 75 grid, computed using the basic MacCormack method with no 
artificial dissipation term. A measure of the width of the mixing layer is the vorticity 
thickness, which is computed based on the mass-weighted velocity profile: 


6 - = 


AU 


d(pu/p) 


dy 

max 


(17) 


where u and pu are spatial averages in the homogeneous x direction. The developed 
structure is shown in contour plots of mixture fraction, density, vorticity and Mach 
number on figures 2a- d. To assess the accuracy of the various schemes a simulation 
with a finer grid was made, hereafter referred to as the high resolution run. Contour 
plots for this case, using a 150 x 150 grid with symmetric TVD scheme (10) and limiter 
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S3, are shown on figure 3. The growth in vorticity thickness is shown on figure 4. 
First the small disturbance grows in the linear regime, and the mean flow grows by 
viscous diffusion. Then the disturbance becomes large and grows non-linearly. It finally 
saturates out when the structure has filled the computational box and the periodicity 
constraint prevents it, or any subharmonics, growing any more. 

The contour plot for the mixture fraction (figure 2a) shows how fluid from the 
upper side (solid contours) and fluid from the lower side (dashed contours) have become 
wrapped around each other. There is a steep gradient formed between the two free- 
stream fluids in the saddle-point region between two structures. The width of this 
strained diffusion layer is the characteristic small length scale of this flow. Contour 
plots of the density and Mach number show that fluid around the top and bottom of 
the vortical structure is accelerated to locally supersonic speeds (peak Mach number 
1.2), and then compresses through nearly- normal shock waves. With the MacCormack 
method we see the development of oscillations, especially on the downstream side of 
the shock wave. There is enough natural viscosity in the flow that the MacCormack 
scheme, with no artificial dissipation term, does not actually fail with the development 
of these weak shock waves. More accurate methods such as spectral methods, would 
show strong oscillations at this point. 

The growth history of the mixing layer computed using the three symmetric TVD 
schemes is shown on figure 4, and compared with the MacCormack computation, and 
with the high-resolution run (150 x 150, S3). It can be seen that there is a wide variation 
between the symmetric limiters. We are in a region of the flow physics where the growth 
rate is very sensitive to the Reynolds number, and any extra damping introduced by 
the numerical method (or ‘clipping’ with a TVD scheme) makes a large difference. 
One could run the simulations at a very high Reynolds number and get less difference 
between the schemes. However in that case the smallest length scales would not be 
resolved by the grid, and the numerical scheme itself would be acting as a sub-grid scale 
model, which defeats the purpose of a direct simulation, where we try to compute till the 
turbulent scales of motion with no modeling errors. We find that all the sy mm etric TVD 
schemes are more dissipative than MacCormack’s method (without an added numerical 
dissipation term), and have lower growth rates and smaller developed structures. The 
least dissipative limiter is S3 followed by S2 and SI. The developed structure for the S3 
limiter is shown on figure 5. The shock wave is well captured in 3 grid points, but the 
developed structure is smaller than with MacCormack’s method. The correct growth 
for this Reynolds number would lead to a structure even bigger than that obtained with 
MacCormack’s method ( c.f . figure 2) 

The upwind schemes on a 75 x 75 grid are compared on figures 6a and 6b. These 
show different characteristics to the symmetric schemes. The start of the non-linear 
growth region is reached earlier in time, which agrees with the higher resolution run. 
However the growth rate in the non-linear region, determined by the slope of the curves 
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at a time of approximately 60, is reduced. In this region the limiter U5 is the least 
dissipative followed by U2, U3 and Ul. In general it was found that the reduction in 
growth rate was tied to the usual descriptor of limiters as compressive or dissipative: 
the more compressive the limiter the sharper the shock waves are captured and, in this 
case, the better the growth rate. The most compressive limiter that was tried was the 
‘superbee’ limiter, U4, which is designed for contact discontinuities. This limiter was 
found to over-predict the growth rate. 

The two methods of implementing the TVD part of the predictor-corrector time 
advance were investigated. Method PCI advances the TVD part at the predictor step 
J7 (2) , while PC2 advances it at step U n (equation (8)). PC2 is TVD for the scalar 
constant coefficient case, but requires more storage than PCI. Plots of growth against 
time are compared on figure 7. At this Mach number there is no significant difference 
between the methods. At lower Mach numbers it was found that PC2 is slightly more 
dissipative than PCI. 

The cost of the methods, relative to the basic MacCormack method, were assessed 
for a 75 x 75 computation run on the Cray X-MP 4/8 at NASA Ames. The symmetric 
TVD schemes were 17% more expensive than the basic MacCormack method for this 
problem, and the upwind schemes were 93% more expensive. The reason for the large 
difference between the upwind and symmetric schemes lies in the ability of the methods 
to be vectorized on the Cray X-MP. For the symmetric schemes with no entropy correc- 
tion all conditional statements can be removed from the TVD loops and these vectorize. 
However the upwind schemes still have to check for a zero in the denominator of equa- 
tion (14c) and this conditional statement prevents vectorization of the main TVD loop. 
Note that an arbitrary small parameter can be added to the denominator of (14c) to 
remove the conditional statement, be we chose to use (14c) to avoid an added parame- 
ter. The simulations of the compressible mixing layer without shock waves presented in 
Sandham (1989) cost about twice as much as the MacCormack method. However those 
simulations had spectral accuracy, and were fully resolved. 

Based on cost and performance for the mixing layer problem the choice of best 
method is between the symmetric TVD scheme with limiter S3, and the upwind schemes 
with limiters U5 or U6. The U6 limiter is a TVB scheme with uniform second-order 
accuracy, and is as good as the best TVD schemes. The upwind schemes U5 and U6 
have the best results for a given grid, but in the current implementation cost 65% more. 

The effect of resolution for the symmetric scheme was evaluated by running four 
simulations with grids of 53 x 53, 75 x 75, 106 x 106 and 150 x 150. The convergence 
of the upwind schemes is expected to be similar. If the scheme were second-order we 
would expect the error to decrease by a factor of 2 for each increase in resolution of y/2. 
The total cost of the simulations increased by a factor of about 3 with each increase in 
resolution. Plots of the growth history for each of these cases are shown on figure 8. 
It can be seen that a large number of points are required if we want to capture all the 
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details of the mixing layer growth correctly. Part of the problem is that of the TVD 
scheme reducing to first-order at points of extrema. The same problem will affect all 
TVD, MUSCL (Van Leer, 1974) and FCT methods which attempt to combine a low 
order dissipative scheme for use around a shock wave, with a high order scheme for 
smooth regions. It is possible that extra accuracy in smooth regions is not worth the 
additional cost since the schemes still reduce in order at extrema, though this issue has 
not been addressed here. Contour plots for the highest resolution case (150 x 150, PCI 
with symmetric TVD and limiter S3) are shown on figure 3. 


5. Discussion 


It is well established that TVD schemes do a good job of capturing shock waves. It 
was the intention of this study to assess the performance of TVD schemes (10) and (14) 
in a flow where other physical processes, besides shock waves, need to be captured. The 
basic finding is that these second-order TVD schemes require very fine grids in order to 
capture all the flow features in the type of physical problem considered here. 

The major drawback with TVD schemes (10) and (14) appears to be the low rate 
of convergence as the mesh is refined, due to a low global order of accuracy. Other 
TVD-type methods which attempt to obtain a higher order of accuracy are high order 
MUSCL (Van Leer, 1974) and FCT (Boris and Book, 1973) schemes. However both 
these methods obtain higher order by increasing the accuracy of the solution only in 
smooth regions of the flow. They still reduce to first-order at points of extrema. Hope 
for future computation of these flows may lie with UNO (Uniformly Non-Oscillatory) 
and ENO (Essentially Non-Oscillatory) schemes (Harten et al., 1987), which try to fit 
high-order polynomials through all parts of the flow, including points of extrema. Very 
impressive results have been obtained with these methods for idealized model problems. 
However the performance of these schemes on more practical flow problems remains 
to be seen. The schemes for reconstruction presented by Harten et al. (1987) do lose 
accuracy at points of extrema, unless accompanied by finer grids. It is not clear what 
the overall gains are, since the operations count and grid stencil for ENO-type schemes 
are larger than for TVD-type schemes. 

Highly accurate numerical methods, such as spectral and Pade schemes, do a very 
good job of resolving turbulent flow features. However when shock waves are in the flow 
a prohibitively low Reynolds number, below that needed for resolution of the turbulence, 
is required in order to resolve the shock waves. Without full resolution these schemes 
rapidly fail. An artificial viscosity term could be added to these methods, provided that 
this term does not affect the computation of the turbulent scales. The dissipation due to 
artificial viscosity, viewed in wave space, would have to be concentrated at wavenumbers 
higher than the largest turbulence wavenumber. 


10 



6. Conclusions 


A class of second-order explicit Total Variation Diminishing schemes have been 
evaluated for numerical simulation of the two-dimensional compressible time-developing 
mixing layer, with the following conclusions: 

(1) The best TVD methods were a symmetric scheme with limiter S3 (equation (12c)) 
and upwind schemes with limiters U5 (equation (15e)) and U6 (equation (16), 
a TVB scheme). The Upwind schemes gave better results but for the current 
implementation were 65% more expensive due to less efficient vectorization on 
the Cray X-MP. Other limiters were more diffusive and it was found that the 
‘compressive’ nature of the limiter was sufficient to predict its performance on the 
mixing layer problem. The less compressive and more dissipative the limiter, the 
slower the growth of the mixing layer. 

(2) All the methods gave slightly different results, and for this problem it is necessary 
to use up to 150 x 150 grid points to properly resolve the growth of a single vortical 
structure. 

(3) At present no ideal numerical method has been found for this flow, which requires 
high accuracy to resolve the growth of an instability in a viscous regime, and the 
ability to handle shock waves without oscillations. There is a need for uniformly 
higher than second-order accurate numerical methods which are efficient and in- 
clude a shock-capturing capability. 
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Figure 1. Schematic of the time-developing mixing layer. 
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Figure 4. Growth versus time, comparing symmetric TVD schemes on a 75 x 75 grid 
with MacCormack method on the same grid, and with a high-resolution run 
on a 150 x 150 grid using the symmetric TVD scheme S3. 
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Figure 6. Growth versus time, upwind TVD schemes on a 75 x 75 grid compared with 
a high-resolution run ( c.f figure 4): (a) Ul, U2, U3 (b) U4, U5, U6. 
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Figure 7. Comparison of predictor-corrector methods PCI and PC2 on a 75 x 75 grid. 
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Figure 8. Effect of increased resolution on growth (TVD S3 method). 
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